/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda - http://marcus.cuda.net
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */
using dnAnalytics.Math;

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal class UserLU : AbstractLU
    {
        public UserLU(Matrix matrix) : base(matrix)
        {
        }

        protected override void DoCompute(Matrix matrix, int[] pivots)
        {
            int order = matrix.Rows;
            for (int i = 0; i < order; i++)
            {
                pivots[i] = i;
            }

            double[] LUcolj = new double[order];

            for (int j = 0; j < order; j++)
            {
                // Make a copy of the j-th column to localize references.
                for (int i = 0; i < order; i++)
                {
                    LUcolj[i] = matrix.ValueAt(i, j); // .Data[indexj + i];
                }

                // Apply previous transformations.
                for (int i = 0; i < order; i++)
                {
                    // Most of the time is spent in the following dot product.
                    int kmax = System.Math.Min(i, j);
                    double s = 0.0;
                    for (int k = 0; k < kmax; k++)
                    {
                        s += matrix.ValueAt(i, k)*LUcolj[k];
                    }

                    LUcolj[i] -= s;
                    matrix.ValueAt(i, j, LUcolj[i]);
                }

                // Find pivot and exchange if necessary.
                int p = j;
                for (int i = j + 1; i < order; i++)
                {
                    if (System.Math.Abs(LUcolj[i]) > System.Math.Abs(LUcolj[p]))
                    {
                        p = i;
                    }
                }
                if (p != j)
                {
                    for (int k = 0; k < order; k++)
                    {
                        double temp = matrix.ValueAt(p, k);
                        matrix.ValueAt(p, k, matrix.ValueAt(j, k));
                        matrix.ValueAt(j, k, temp);
                    }
                    pivots[j] = p;
                }

                // Compute multipliers.
                if (j < order & !FloatingPointUtilities.EqualsWithinSignificantFigures(matrix.ValueAt(j, j), 0.0, 15))
                {
                    for (int i = j + 1; i < order; i++)
                    {
                        matrix.ValueAt(i, j, (matrix.ValueAt(i, j)/matrix.ValueAt(j, j)));
                    }
                }
            }
        }

        protected override Matrix ComputeInverse(Matrix matrix, int[] pivots)
        {
            int order = matrix.Rows;
            Matrix inverse = matrix.CreateMatrix(order, order);
            for (int i = 0; i < order; i++)
            {
                inverse.ValueAt(i, i, 1.0);
            }

            Solve(mMatrix, mPivots, inverse);
            return inverse;
        }

        protected override void Solve(Matrix factor, int[] pivots, Matrix result)
        {
            Pivot(result);
            int order = factor.Rows;
            // Solve L*Y = B(piv,:)
            for (int k = 0; k < order; k++)
            {
                //                int korder = k * order;
                for (int i = k + 1; i < order; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        double temp = result.ValueAt(k, j)*factor.ValueAt(i, k);
                        result.ValueAt(i, j, result.ValueAt(i, j) - temp);
                    }
                }
            }
            // Solve U*X = Y;
            for (int k = order - 1; k >= 0; k--)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    result.ValueAt(k, j, (result.ValueAt(k, j)/factor.ValueAt(k, k)));
                }
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < result.Columns; j++)
                    {
                        double temp = result.ValueAt(k, j)*factor.ValueAt(i, k);

                        result.ValueAt(i, j, result.ValueAt(i, j) - temp);
                    }
                }
            }
        }

        protected override void Solve(Matrix factor, int[] pivots, Vector result)
        {
            Pivot(result);
            int order = factor.Rows;
            // Solve L*Y = B(piv,:)
            for (int k = 0; k < order; k++)
            {
                for (int i = k + 1; i < order; i++)
                {
                    result[i] -= result[k]*factor.ValueAt(i, k);
                }
            }
            // Solve U*X = Y;
            for (int k = order - 1; k >= 0; k--)
            {
                result[k] /= factor.ValueAt(k, k);
                for (int i = 0; i < k; i++)
                {
                    result[i] -= result[k]*factor.ValueAt(i, k);
                }
            }
        }

        private void Pivot(Vector data)
        {
            for (int i = 0; i < mPivots.Length; i++)
            {
                if (mPivots[i] != i)
                {
                    int p = mPivots[i];
                    double temp = data[p];
                    data[p] = data[i];
                    data[i] = temp;
                }
            }
        }

        private void Pivot(Matrix data)
        {
            for (int i = 0; i < mPivots.Length; i++)
            {
                if (mPivots[i] != i)
                {
                    int p = mPivots[i];
                    for (int j = 0; j < data.Columns; j++)
                    {
                        double temp = data.ValueAt(p, j);
                        data.ValueAt(p, j, data.ValueAt(i, j));
                        data.ValueAt(i, j, temp);
                    }
                }
            }
        }
    }
}